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COMPARISON OF NUMERICAL TECHNIQUES FOR THE EVALUATION OF 
THE DOPPLER BROADENING FUNCTIONS ty(x, 8) AND xUG) 
by R. Bruce Can right, Jr., and Thor T. Semler 
Lewis Research Center 

SUMMARY 

The two Doppler broadening functions, ^/(x, Q) and x( x > 0)> are necessary for the 
computation of accurate temperature dependent resonance neutron cross sections. With 
the present use of many resonances and ever increasing resonance energies for reactor 
physics calculations in support of fast reactor development, the direct numerical quad- 
rature evaluation of these functions has become prohibitively expensive in terms of com- 
puter time. 

In this report direct numerical quadrature is compared with the techniques of Gaus- 
sian quadrature, contour integrations, and cubic spline interpolation in terms of both 
accuracy and ease of computation. A form of contour integration due to A. M. Turing 
(1943) is found to be much faster than (approximately six times) and as accurate as 
direct numerical quadrature. 

Execution times for the evaluation of these functions on the IBM 7094-11 computer 
and FORTRAN IV listings are included. 


INTRODUCTION 

It was noted early in the study of neutron physics that the neutron cross section, or 
interaction probability, of materials is dependent on the temperature of the material. 
Thus, a rapid variation of cross section with respect to neutron energy could be inferred 
(refs. 1 and 2). It was then pointed out by Breit and Wigner (ref. 3) that the nuclear 
cross sections for the formation of the compound nucleus, or target nucleus plus neu- 
tron, should exhibit a resonance structure. For the reaction, neutron in and gamma 
out, the cross section <x should vary as a function of neutron energy E as shown in 
the following equation (ref. 4): 



( 1 ) 



where k is the wave number of the neutron in the center of the mass system, gj is the 
spin factor (2J + 1)/2(2I + 1), J is the spin quantum number of the compound system and 
I that of the nucleus, T n is the scattering width, is the radiative width, E Q is the 
energy of resonance, and T is the total width. Equation (1) is correct for stationary 
nuclei. However, since the nuclei of all materials are in thermal agitation, equation (1) 
must be modified to account for the effect of thermal motion. Briefly, the nuclei are 
assumed to have a Maxwellian velocity distribution and the cross section is folded into 
this velocity distribution. The resulting cross section is then a function of temperature 
as well as neutron energy; this cross section is the Doppler broadened cross section. 
Two functions arise from the derivation of the Doppler broadening of the cross section; 
they are i//(x, 0) and \(x, 0). The derivation of the function j//(x, 0) is shown in detail 
in appendix A. The respective definitions are as follows: 



Here x = (2/r)(E - E Q ), 0 = r(4mTE o /M)~*/ 2 , m is the neutron mass, M is the nuclear 
mass, and T is the effective temperature in eV T s. 

Plots of i//(x, 0) and x(x, 0) are shown for selected values of 0 and a range of x 
from 0 to 10 in figures 1 and 2. Detailed plots of the regions of interest are shown in 
figures 3 and 4. 

In the past various forms of ’’brute force’’ direct numerical quadrature have been 
used to evaluate \p and \ (ref. 5). However, with the present use of many resonances 
and ever increasing resonance energies for reactor physics calculations in support of 
fast reactor development direct numerical quadrature has become prohibitively expen- 
sive in terms of computer time. The following report evaluates various techniques 
which might avoid the excessively long computer runs needed to calculate the Doppler 
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Figure 3. - Portion of X, 8) surface. Figure 4. - Portion of (x X, 8) surface. 


broadened cross sections by direct numerical quadrature. Such an evaluation is valuable 
because of present widespread interest in such fast reactors, and because "lengthy and 
involved routines giving wrong values for the (Doppler broadening) functions are preva- 
lent" (ref. 6). 

It should be noted that, although these integrals have been extensively studied (refs. 

6 to 12), there is virtually no standard notation for either the Doppler broadening func- 
tions or the independent variables. We follow the notation of Dresner (ref. 4). Here 9 
is proportional to T” 1//2 and x is proportional to the energy distance from the reso- 
nance peak. Typical transformations for comparing our notation to others include 


Dresner Bhat and Lee- Whiting Matta and Reichel 



NUMERICAL TECHNIQUES 

In order to evaluate temperature dependent continuous neutron cross sections, one 
must evaluate and x many times; therefore, both accuracy and as much speed in 
computation as possible are required. Special series methods, numerical quadrature, 
contour integration techniques, and interpolation have been considered for their evalua- 
tion; these studies indicate that a specialized quadrature and contour integration yield 
the accuracy required. 


v 


Asymptotic Series and Limiting Cases 


Several special forms for and \ are known which are simple to evaluate, but 
taken all together, do not cover much of the physically meaningful range of x and 9. 
For most work |x| < 1000, 0.001 < 6 < 100. For this reason such special forms are 
used only as boundary check values for the other more general methods. 

For example, for x9 >> 1 

if/ (x, 6) = — - — (1 + — - 3 - 2 -" -— + . . . ) (ref. 4) 
i + x 2 \ e 2 (i + x 2 )2 / 


(ref. 4) 


(ref. 4) 


Direct Numerical Quadrature 

For large-scale problems, several hundred resonances and several thousand energy 
points, the standard quadrature techniques in use at Lewis Research Center are too slow 
by at least two orders of magnitude. Because of the form of the integrand, programs 
using a Gaussian quadrature and also a Simpson’s rule quadrature containing step modi- 
fication logic were tried. This adaptive Simpson's rule quadrature, called method A, 
is coded in FORTRAN in appendix B as function SIMPS 1. Both programs produced large 
errors for some values of x and 9; the critical parameter is the product x0. Both the 
shape of the integrands and the choice of limits of integration depend on this product. 

Let us transform the integrands and examine their behavior. Let 


t = - (y - x) (6) 

2 


Also, 


^ (x, °° ) = • 


X(x, °° ) = ■ 


1 + x^ 


2x 

l + x 2 


\p(0,e) =-i x(0, e ) 


Then 


5 



( 7 ) 



-t 

The damping factor e in these integrands is now independent of x and 9 both 
(with a maximum at t = 0) and the remaining terms are 


f s 



(9) 



( 10 ) 


These terms are sketched in figures 5 and 6. Note that both f and g have critical 
points which depend on -x9/2. 
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The integrals for xj/ and 
Hermite quadrature (ref. 13). 


X are now in convenient form for application of Gauss- 
That is, 


* 


X = 



2 N 
f(t) dt ^A^t) 

i=l 

N 

g(t) dt ^A-gfr.) 

U1 


(ID 


( 12 ) . 


where the t^ and A^ are tabulated for many N; for example, we chose N = 16 and 
N = 32. 

However, it can be seen from figures 5 and 6 that, for Gauss- Hermite quadrature, 
no matter what t. are picked, there will be both x and 9 such that much of the area 
under f and g is far from any of the t. . For example, for N = 16, the maximum 
b = ±4. 689 - meaning f and g are sampled only in this range; typical values of x = 
20, 6-2 place the critical points at t = -20. This inability to sample from significant 
values of f(t) and g(t) explains why Gauss- Hermite quadrature fails to give accurate 
results. 

On examining figure 6 from the point of view of using a Simpson’s rule quadrature, 
it can be seen that, for some values of x and 9, the x integrand has a positive por- 
tion and a negative portion. To achieve and maintain accuracy, . x must be broken into 
positive and negative portions and evaluated separately. With these modifications, the 
adaptive Simpson's rule was successfully applied to \fs and xi this shall be referred to 
as method A. One should note that the two parts -<*> < t < -xB/2 and -xd/2 < t < +°° 
of x sometimes agree in magnitude to 12 significant figures. For comparison, another 
quadrature method in use at this center called method D was also studied. 


Contou r Integration 


Integrals of the form 




dt, as we see in equations (7) and (8) can 


sometimes be evaluated by contour integration in the complex plane (ref. 14). A. M. 
Turing first suggested this method for evaluating \(/ and x m 1943 (ref. 15). The 
Turing method is attractive because it should require much less arithmetic than numer- 
ical quadrature. This method has been formulated independently by Bhat and Lee- 
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Whiting, and by Matta and Reichel. We have applied the method of Bhat and Lee- Whiting 
(ref. 6), called method B and the similar method of Matta and Reichel (ref. 12), which 
we call method C. These methods are very useful for many related integrals as well 
(ref. 12). Methods B and C use series with a parameter h (and corresponding error 
estimate E(h)) which can be varied to achieve desired accuracy. In addition to the ser- 
ies contribution, both methods have a contribution from the poles (places where the de- 
nominator of the integrands is zero) which depends on whether 6 is inside, on, or out- 
side the contour chosen for integration. The formulas are the following: For method B 
let a s x0/2, b s q/ 2, and recall transformations (4) and (5). Then 


^ (a, b) 


hb 

“00 

y e-" 2 '' 2 

9 9 cos (2ab) - e 2?rb / h cos f 2ffa - 2ab\ 

+ 2P(h)e^ b _a ) ' b ' 

v; 

/ J (a - nh) 2 + b 2 

D 


- “ 00 



+ real part of E(h) (13) 


X (a,b)=-^ 


frr 


" oo 

I 


2.2 

-n h 


(a - nh) 


(a - nh) 2 + b 2 


L - co 


4P(h)e 


9 9 sin (2ab) + e 2nb/h sin l 2 ^ - 2ab 

(b -a 2 ) \ h , 


D 


+ imaginary part of E(h) (14) 


where 


D = 1 - 2e 2?rb/h cos^ + e 4 ^/ 11 

h 


and 


P = 0 for 2 77/h < G 


P = 1/2 for 27r/h -dr (ref. 6) 


P = 1 for 277/h > 6 J 


Finally, 
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E (h) 1 


Va 2 + b 2 e~ n2/h 



2 9 
1_ - b 2 



h 2 



(15) 


This method is coded in FORTRAN in appendix C, as subroutine PSCH. 

o 

Similarly, for method C, let t = 1/8 . Then 


• / . \ h 2h 

(x, t) = + 


V 2 2 

e~ n h (1 + x 2 + 4tn 2 h 2 ) + pp< _ E(h) (16) 


V; ( l + x 2 ) yTn^(i_ x 2 + 4tn 2 h 2j + 4x 2 * * 


X(x,t) =■ 


where 


2hx 

4hx \ 

e" n h (1 

+ x 2 - 

4tn 2 h 2 ) 0 

(1 + x 2 ) 

v ; i 

n: 

^( 1 - x2 + 

4tn 2 h 2 

/ ♦ 4x^ 

o 

II 

for 

27r/h < 9 

or 

t < h 2 /4;r 2 

P - 1/2 

for 

2n/h = 8 

or 

t = h 2 / 

P = 1 

for 

2 n/h > 8 

or 

t > h 2 /4?r‘ 


-2jpQ 1+ 5E(h)J (17) 


Also, 


and 


P t = a /AC_-_BD\ = g BC ± AD) 

k C 2 + D 2 / \C 2 + D 2 


a _ jl e - (x 2 / 4t+7r/hVt - 1/ 4t) 
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A = cos — 
2t 


B = sin — 
2t 


C 


= e - ff /hVt _ cosJ2L 
hVt 


D = sin 


JTX 

hVt 


Finally, 


E(h) | < 


2 ,. 2 

r TT /h 


«(l - .-*/»* 


(18) 


This method is coded in FORTRAN in appendix C, as subroutine PSIPHI. 

Using these estimates for the error, and knowing the accuracy desired, one can 
choose the largest possible h which will achieve this accuracy. For method B, h = 1. 0 
was used; for method C, h = 0. 75. The series of each method was truncated at nine 
terms. The same number of terms is used in order to make comparisons of execution 
times. 


Interpolation 

Using methods A, B, and C tables were built of and x for various x and 6 in 
hopes that interpolation would be possible. These tables were fit with cubic splines 
(refs. 16 and 17). Briefly, cubic splines are cubic polynomials passed piecewise through 
each pair of adjacent data points. The coefficients of these cubics are chosen so as to 
match first and second derivatives at the data points; this is usually called the spline 
property. These points were used to generate contour plots as shown in figures 7 and 8. 
This method of interpolation has not yet achieved the accuracy of methods A, B, or C. 
Interpolation remains attractive and is a subject for future work. 
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RESULTS 


Accuracy Comparisons 


The Gaussian quadrature techniques as used herein do not give satisfactory relative 
accuracy in computing i// and x- This is true whether or not the weighting function 
_t^ 

e is removed from the integrand. For most values of x and 9 this method is satis- 
factory, but wrong results occur in an unpredictable fashion. Since the error estimate 
associated with Gaussian quadrature assumes the use of functionals not readily obtain- 
able, it has not been used. The results are in error because this quadrature uses fixed 
abscissas, while, as discussed in the section Direct Numerical Quadrature, the peaks of 
the integrands for \p and x vary as the product x9. 

The adaptive Simpson's rule method (method A), with numerical results given in 
table I, yields satisfactory accuracy for our work, four significant figures. This meth- 
od was the first one found to be sufficiently accurate; it was used as a check against 
published values and against other methods discussed in this report. 

The contour integration methods, methods B and C, have accuracies controllable 
by choice of a parameter h. See tables n and III. They both have corresponding error 
estimates E(h). For a given fixed h, methods B and C give much the same results, 
and speed becomes the deciding factor. However, there is one region where method B 

TABLE I. - TYPICAL VALUES OBTAINED BY METHOD A 


(a) Values of i // 


X 

e 


0. 01 

0. 5 

1 

2 

0. 001 

0. 00881246 

0.341353 

0. 545648 

0.757895 

. l 

. 00881246 

.341191 

. 544852 

. 755166 

l 

. 00881224 

.325577 

. 472498 

. 540149 

4 

. 00880896 

. 164687 


. 0643058 


(b) Values of x 


X 

8 

0.01 

0. 5 

1 

2 

0. 001 
. 1 
I 
4 

0.991191X10' 7 
• 991191X10' 5 
. 991175x10" 4 
.396371X10" 3 

0. 164663X10" 3 
. 0164612 
. 159261 
. 400753 

0. 454362x10" 3 
. 0453743 
. 408498 
. 504762 

0. 968510X10" 3 
.0965705 
.738117 
. 481550 
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TABLE II. - TYPICAL VALUES OBTAINED BY METHOD B 


(a) Values of \p 


X 

e 


0 . 01 

0.5 

1 

2 

0.001 

0 . 00881269 

0. 341350 

0.545639 

0.757862 

. 1 

. 00881239 

.341189 

. 544844 

. 755133 

1 

. 00881233 

. 325574 

. 472496 

.540137 

4 

. 00880900 

. 164686 

. 0915985 

. 0643081 


(b) Values of x 


X 

e 

0 . 001 

0 . 5 

1 

2 

0 . 001 
. 1 

1 

4 

0 . 995718 x 10 " 7 
. 988693 x 10' 5 
. 992042 x 10' 4 
. 396523 X 10' 3 

0 . 164663 x 10' 3 
. 0164608 
. 159267 
.400761 

0 . 45 43 63 x 1 O ' 3 
.0453872 
. 408535 
. 504784 

0 . 968537 x 10' 3 
. 0965732 
.738138 
.481569 


TABLE III. - TYPICAL VALUES OBTAINED BY METHOD C 
(a) Values of 


X 

e 


0 . 01 

0 . 5 

1 

2 

0.001 

0. 00881246 

0.341351 

0.545641 

0 . 757872 

. 1 

.00881249 

.341189 

.544846 

.755143 

1 

. 00881225 

.325575 

. 472498 

. 540145 

4 

. 00880898 

. 164687 

. 0915994 

. 0643072 


(b) Values of \ 


X 

e 

0.01 

0.5 

1 

2 

0 . 001 
. 1 

1 

4 

0 . 990694 X 10' 7 
. 991765 X 10' 5 
. 990592 x 10' 4 
. 396352 X 10' 3 

0 . 164662 x 10' 3 
. 0164607 
. 159265 
. 400756 

0 . 454359 X 10' 3 
.0453867 
. 408530 
. 504770 

0 . 968511 X 10" 3 
. 0965709 
.738117 
. 481536 









































TABLE IV. - TYPICAL VALUES OBTAINED BY METHOD D 


(a) Values of i// 


X 

e 


0. 01 

0.5 

1 

2 

0. 001 

0. 00881212 

0. 341334 

0.545639 

0.757872 

. l 

. 00881228 

.341179 

. 544870 

. 755142 

l 

. 00881233 

.325574 

. 472496 

.540137 

4 

. 00880896 

. 164684 

. 0916000 

.0643058 


(b) Values of x 


X 

e 


0. 01 

0.5 

1 

2 

0.001 

0.951862x10' 7 

0. 164369x10" 3 

0. 455154x10" 3 

0. 968535X10' 3 

. 1 

. 941018x10' ^ 

.0164325 

. 0454529 

.0965704 

1 

.994859x10' 4 

. 159274 

. 408537 

. 738117 

4 

.395791X10' 3 

. 400745 

. 504769 

. 481550 


loses all accuracy, while method C remains stable. For 6 ~ 0. 001 (high energies 
> 1 eV) method B failed, even with h = 0. 5 and 20 terms in the series. In this region 
method C, with h = 0. 75 and nine terms as before, gives four significant figures for 
i// and one for We conclude that, for small 9, method B should be used with caution. 

Some sample results obtained by methods A, B, C, and D are given in tables I to 
IV. 


Speed Comparisons 

Because the numerical quadrature of method A requires far more arithmetic oper- 
ations than methods B and C, method A has been modified to split the argument ranges 
into two regions and to use an asymptotic expansion when possible. Values for x were 
chosen uniformly in the range 0 < x < 1000. (Note that i^(-x) = i//(x) and x(-x) = - \(x ). ) 
Two sets of 6 ranges were used: in the range 0 < 9 < 1, the pole contributions are 
always present; while in the range 0 ~ 9 < 100, these contributions are usually zero. 
Execution times on an IBM 7094 for 2000 evaluations are given in table V. An example 
of a Doppler broadening resonance cross section code which may be used with the sub- 
routines herein is shown in appendix D. 
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TABLE V. - COMPARISON OF EXECUTION TIMES 


8 range 

Method 

A (integrals only) 

A (integrals + series) 

B 

C 

D 

Time, sec 

0 < 9 < 100 

a 5. 30 

2. 17 

0.77 

1.87 

2. 16 

0 < 8 < I 

— 

11.5 

2. 05 

2.20 

— 


a Run for comparison only--it was known that this method would not be 
fast enough to be practicable. 


Discussion of Execution Times 

Several conclusions can be drawn from table I. First, even the slow numerical 
quadrature method is accurate, and as mentioned earlier "lengthy and involved routines 
giving wrong values for the i p and \ functions are prevalent" (ref. 6). Study of these 
integrals began because an existing routine was very slow, about six times slower than 
our quadrature. Note that adding a test and an asymptotic expansion to method A in- 
creases the speed nearly 2^ times. Most important for this study, method B is nearly 
2^ times faster than method C for the distribution 0 < 8 < 100. This is because terms 
in the series of method B can be written more compactly than those of method C; and 
only for small 8 is the pole contribution needed; that is, the series is the only calcu- 
lation required. Fourth, for 0 < 8 < 1 most of this speed advantage is gone. Now the 
pole calculations are always needed. Method B, coded as subroutine PSCH listed in ap- 
pendix C, is still to be preferred, for methods B and C are about equally accurate. 


CONCLUDING REMARKS 

Direct numerical quadrature using an adaptive Simpson's rule technique yields re- 
sults accurate to four significant figures for both ^/(x, 8) and \( x > #)• However, this 
method is quite slow and expensive in terms of computer time. 

The method of Gaussian quadrature used herein was evaluated and found to be er- 
ratic in terms of error and somewhat slower than the techniques of contour integration. 
Error estimation for the Gaussian quadrature proved to be incapable of rapid calculation. 

Both contour integration methods, method B using subroutine PSCH, and method C 
using subroutine PSIPHI gave very rapid and accurate results. The accuracy is com- 
parable with that of direct numerical quadrature. 
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The method of cubic spline interpolation failed to achieve the accuracy of the other 
methods evaluated. 

Of the numerical methods considered and evaluated contour integration, subroutines 
PSCH and PSIPHI, achieved the greatest speed of computation combined with quite satis- 
factory accuracy. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, April 10, 1972, 

132-80. 



APPENDIX A 


DERIVATION OF vy(x, 0) - THE DOPPLER BROADENING FUNCTION 

If a beam of neutrons with velocity v impinges upon a group of nuclei at rest, the 
reaction rate for process x is proportional to vcr (v), where a (v) is the nuclear cross 
section for process x as a function of velocity. However, the nuclei are not stationary 
but are in thermal motion. We may assume the nuclei have a Maxwellian velocity dis- 
tribution characterized by a temperature T; then the three-dimensional Maxwell veloc- 
ity distribution may be written as 


P(V) dV 


=(— ] 

\27rkT/ 


3/ ' 2 e -MV 2 /2kT dy dy dy 


where T is the effective temperature, M is the nuclear mass, and V is the nuclear 
velocity. 

The value of the velocity of the neutron, of mass m, relative to the nucleus is then 
the magnitude of the difference of the velocity vectors. That is, 


v 


r 


= V - v 


12 — — 

The corresponding energy is then E r = g mv r> Since P(v ) dv is the probability of a 
nucleus having a velocity within dv about v , the probable number of reactions to type 
x per second is 


v r a x( E r) P ^ 


The total probability is then 


a x (E) =i f v r a x (E r )P(V) dV 

If a coordinate system is chosen such that the z- direction is parallel to the neutron 
velocity, then 


E =— m(v - V) 2 
r 2 


1 

2 


m 


(v- 


V 
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Since the neutron velocity is quite high compared to the thermal motion of the nuclei, 
2 2 2 

the V_, V . and V, terms may be neglected. Hence, 
z. x y 

E = - m(v 2 - 2vV ) 
r 2 z 


Thus 


v r 


Integrating out the V x and V y velocity components from the Maxwellian results 
on the following distribution for the z- component: 




P(V) dV x dV y 


/ M \V 2 e -MV z /2kT 

\2ffkT/ 


The Breit-Wigner expression for a single- level resonance cross section for absorp- 
tion is 


a (E) = a 
y v ' o 




4(E - e 0 ) 2 + r 2 


o 

where a Q = 477 (r n /r), h is neutron wave length over 277, g is the statistical 

weighting factor, T n is the scattering width, T is the total width, T y is the radiative 
width, and E Q is the resonance energy. Substituting into the previous equation gives 



If one defines 


x S l(E-E 0 ) 
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y'p< E r- E o) 


A / 4mkTE , r^o m . x 


(a very good approximation for large E q ) and 


0 = r 


- o 


— J--- 

r I/ e 


exp - 0 2 (y - x) 2 ] dy 

1 + y 2 


where 


a (E) = a — il/(x,9) 


r i e 


V/(x, 0) = ■ 


/ 


exp - - 0 2 (y - x) 2 dy 
. 4 J 

1 + y 2 


The function 


X(x, 0) = ■ 


/ 


2y exp - - 0 2 (y - x) 2 dy 
L 4 . 

1 + y 2 


may be derived in the same manner using the same assumptions. 
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APPENDIX B 


FORTRAN IV LISTING OF SIMPS1 


This function SIMPS 1 is a Simpson's rule integrating routine, which integrates the 
external function FUNC1 from XMIN to XMAX. When sufficient accuracy is not attained, 
the flag KER is incremented by 1. 


FUNCTION SINPS1 ( XI* IN, XMAX, FUNCl.KER ) 

DIMENSION V I 200 ) »H I 20 0 ) » A ( 200 ) , B I 2CG ) » C l 2QC ) * P I 2CC ) » E I 2C0 ) 
DATA T/3.CE-4/ 

IFIXUN.EC.XPAX ) GO TO IB 
V ( 1 ) =XM I N 

H I 1 ) =0 . 5 < (XKAX-Xf'IN) 

A ( 1 ) = FUNC 1 1 XN 1 IN ) 

8 ( 1 ) = F U N C 1 ( XWIN + F ( 1 ) ) 

C ( 1 )=FUNC 1 ( XMX ) 

P ( 1 )=H< 1 >* ( A ( 1) +4 ,C*E (1) +C ( 1 ) ) 

EI1)=FI1) 

A N S = P < 1 » 

N= 1 

FRAC=2.0*T 

1 FRAC=C.5«FRAC 

2 TEST = ABS(FPAC*ANS ) 

K = N 

2 CC 7 1=1, K 

4 IF (AES( E ( I ) ) .LE.TEST) GO TO 7 

5 N = N + 1 

VI N ) = V I Il+MI) 

H I N ) =C • 5 < H I I ) 

A I N ) = E I I ) 

BIN ) = FUNC 1 ( V ( N) +MN) ) 

CIN )=C( I ) 

PIN ) =H(N)* ( A ( N ) +4.0*e IN ) +C I N ) ) 

G= F I I ) 

H I I I = H ( N ) 

ei I ) = FUNC 1 (VI !)+(-( I ) ) 

C I I ) = A ( N ) 

PI I ) =H I I ) * I A ( I)+4.C*ei I ) +C I I ) ) 

G= P I I l + PIN )-Q 
ANS=ANS+C 

E I N ) = C 

£ IF IN-200) 7,13, 13 
7 CONTINUE 
£ IF (N-K) 9,9,2 
9 Q = C.O 
1C CC 11 1=1, N 

11 Q = G + E 1 1 ) 

12 IF I AES (C )-T* ABS I ANS ) ) 14,14,1 
13 KER=KER+1 

1 A A N S = C . J 

15 CC 16 1 = 1, N 

16 ANS=ANS + F I I ) 
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APPENDIX C 


FORTRAN IV LISTINGS OF METHODS B AND C 

Subroutine PSCH uses the technique of Bhat and Lee- Whiting (ref. 6). The argu- 
ments are XX and TT where 

XX natural distance X 
TT temperature parameter 

The values U and V are returned where 

U first Doppler broadening integral 
V second Doppler broadening integral 


SlBFTC PSCHEL 

C FPCM EH/ST AND LEE-WHITING, P.278 

C H=l. 

SLERCUTINE PSCH ( XX , TT , U , V ) 

DIPENSION EPN(5 ) ,EN2( 10) 

DATA FI22/12. 566371/ 

DATA RTPI/1. 7724539/ , P IH/ 1 . 5707963/ 

DATA PI , F I 2 , N /3. 1415927, 6. 2831854, 5/ 

DATA (EPM I ) , 1 = 2,5) / . 1 1 709966E0 , . 583C04 89E-2 , . 392E 2 561E-4 , 

1 .35821C59E-7 / 

DATA (EN2 U ) , 1=2,10) /.24785S99, .11709966, .33549615E-1 ».58300489E- 
*2,.61448264E-3, . 39 282 56 1 E-4, . 1 5 2 3 1 502E-5 , . 3 582 1059E-7 , . 5 1095996E-9 
*/ 

C 

X C = A E S ( XX) 

THETA = AES (TT ) 

Y=THETA/2 . 

X=XC*Y 

C ARE WE IN DANGER 

XTEST=X-A INT ( X) 

IF ( (XTEST.LT. .01 .OR. XTEST.C-T..59) .AND. A BS ( Y ) . L T . . 01 ) GO TO 40 
IF (XC.GT ,75.*SQRT ( l.+l. / ( Y*Y ) ) ) GO TO 30 
IF(Y.LT..C01 ) GO TO 40 

C - BEGIN SERIES SUPATION -- - - 

Y 2 = Y * Y 

C SET N = G TERM — NOTE PI IS IN THIS ONE 

D = P I* (X + X + Y2 ) 

SLi = Y / C 
SV=X/C 
DC 1C 1=2, N 
AP = I - 1 
XNF=X-AP 
X N P = X 4 A P 
DP=XNP*XNP4Y2 
DP=XNP«XNP4Y2 



SU = SU 4 EKNII ) *Y* ( 1«/0P+1./DM) 

SV = SV 4 EVN(I)*(XNP/DP4-XN! y l/DM) 

1C CONTINUE 

IF(Y.GT.FI) GO TO 20 
P=2 • 

IF(Y.EQ.FI) P=l. 

C EEGIN POLE CCNTR I RUT ION 

XP2=X*P 1 2 
EYF2 = EXP ( Y*P I 2 ) 

15 CONTINUE 
XY2=X*Y*2. 

SXY2=SIN(XY2) 

CXY2=COSIXY2) 

SXP2=SINIXP2) 

CXP2=C0S ( XP2 ) 

EYX=EXP (Y2-X*X) 

0=1.- EYP2*( 2.*CXP2-EYP2) 

St = SU 4 P*EYX*(CXY2-EYP2*(CXP2*CXY2 + SXP2*SXY2 ) ) /C 
SV = SV - P*EYX*(SXY24-EYP2*( SXP 2*CXY2-CXP2* SXY2 ) ) 'C 
2C U=SU*Y*RTPI 

V=SV*THETA*RTPI 
IFIXX.LT. 0.) v=-v 
GC TO 5C 
3C F = PT P I * Y 

SL= 1 • / < (l. + XC*XC I*F) 

SV=SL*XC 
GC TO 2C 
AC CONTINUE 

IFIXC.LT. A. .AND. THETA. LE..C09) WRITE(6,1C1) XD, THETA 
1C1 FCRKAT I AH eAD,2GlA.6> 

C THIS IS EXPANSION WITH H=C.5 

Y2=Y*Y 

C=PI*(X*>+Y2) 

SL=Y/C 

SV=X/C 

DC A 5 1=2,10 
AH=0 • 5*Fl 0 AT ( 1-1 ) 

XNP= X-AK 
X N H = X 4 A H 
CP = XNP*XH P4Y2 
CP = XNK*XN P4Y2 

SU = SU 4 EN2II )*Y*(1./OP+1./OM) 

SV = SV 4 EN2 ( I )*(XNP/OP+XNH/DM) 

A 5 CONTINUE 
SU=SU*0.5 
SV=SV*0 .5 

IF(Y.GT.FIH) GO TC 2C 
P = 2 . 

IF(Y.EQ.FIH) P=l. 

XP2 = X*P 12 2 
EYP2=EXP(Y*PI22) 

GC TC 15 
5C RETURN 
ENC 


24 



oono o non oooo 


Subroutine PSIPHI uses the technique of Matta and Reichel (ref. 12). The arguments 
and values are the same as for routine PSCH. 


SIBFTC NATTA 


FRCM NATTA+RE ICHEL MATH. COMP. P.34C APRIL 1971 
H= 1. , 11 TERMS (FOR COMPARISON) 

SUE R CUT IN‘E PS IPH I ( XX , TT ,U, V ) 

OATA H/C.75/ 

DATA PI/2.14159265/ , RTP I / 1 . 7724 539 / 

DIMENSION ENH (9 ) , FNH ( 9 ) 

DATA FNH/ 2. 25,9. ,20.25,36. ,56.25,81. ,110.25,144. ,182.25/ 
DATA ENH / . 5697 8 2 83 , . 105399 2 2 , . 632 97 1 54E-2 , . 12 34C5 81E-3 , 

* . 78 11485 4E-6,. 16052281E-8, . 107C 92 3 2E- 1 1 , . 23 155 228E- 1 5 , 

* .16310129E-19 / 

DATA BOUND/. 14248292E-1/ 

X = ABS (XX ) 

TH = AES ( TT ) 

T=1.C/TH<*2 

X2=X*X 

START SERIES — SET N=0 TERM 

XP=1.0+X2 
XN=1.C-X2 
SU s C . 5/ XF 
SV=SU 

X24=4,0*>2 

DC 1C N= 1 , 9 

FNHT=FNH(N)+T 

DENCM=X24 + ( XN+FNHT ) **2 

TERM=ENH(N)/DENOM 

SU = SU + TERM* ( XP+FNHT ) 

SV = SV + T ERM* ( XP-FNHT ) 

1C CONTINUE 


SU=SU*2.C*H/RTPI 

SV=SV*2.C*H/RTPI*X 

START POLE TESTS 

IF (T .IT .ECUND) GOTO 20 _ . . 

RTT=SCRT(T> 

EXPCN=XN*Q.25/T - PI/H/RTT 
IFlEXPON.LT. (-34.)) GO TO 3C 
ARG1=X*0.5/T 
ARG2=PI*X/(H*RTT ) 

A=CCS ( ARC 1 ) 

B=S I N ( ARC 1 ) 

C = EXP(-PI/(H*RTT) ) - COS ( ARG 2 ) 

D*S I N ( ARC 2 ) 

TERM=RTPI/RTT*EXP(EXPON)/(C*C+D*D) 
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APPENDIX D 


A SAMPLE DOPPLER BROADENING CROSS SECTION CODE 
WHICH REQUIRES VALUES OF AND x 


$IBF7C RP3TAP 
C 

LCC-ICAL START 
DIMENSION CATUM( 12003) 

Cl REN SI ON CATUMS ( ^001 ,3 ) 

EQUIVALENCE ( DATU V ( 1 ) , 0 ATUMS ( 1 , 1 ) ) 

DIRENSICN T ITLE ( 24 ) » UGI150), DUI15C) 

DIRE N SIGN GG ( 20 ) 

Cl REN SIGN EC ( 300 ), GARNI 300 ) ,GAMG( 300 ) , GAR ( 300 ) , GARE ( 300 ) 

DIR ENSIGN SIGCOI 300) *SIGF0(3C0) ,SIGSOt 300) ,SIGT0(3CC) ,C0N3(300) 
0 IRENS I CN A(300),SIG1(300),SIG2(3C0),SIG3(300) 

Cl RE NS I ON X (300) , THETA ( 300 ) , PSII3C0) , PHII3C0) 

C 

C 

CCRRCN NPES, X, THETA, PSI , PHI 


START*. TFUE • 

N U N IT = 1 
REWINC NLNIT 
CALL SKF ILE ( NUN I T ) 

BACKSPACE NLNIT 
1 CONTINUE 

I F ( ST ART ) GC TO 111 
C FIND ENC OF THIS TAPE 

C PUT CN IC RECORD 

WRITE (NUN IT ) ID 
C WRITE GCCC STUFF 

JWCRCS=3*NWCR0S 
WRITE (NUNIT )DATUK 

75 FCRRAT( A6 ,4X, 110 > 

76 FCRR AT ( IF 1 » A6 ) 

77 FCRRAT ( 1CX , 3 F 13 .4 ) 

111 CONTINUE 

START*. FALSE. 

RE AC (5,50) ( TITLE ( M ) ,M= 1 , 24 ) 

RE AC (5,75 ) IC.NWCRCS 

READ (5,51) T,ANU,UKT,SIGCT,SIGFT,USTART,NRES,NDU 
RE AC (5,53) SIGMAT,SIGPOT,G 
READ (5,55) ( UG ( K ) , DU ( K ) , K= 1 , NDU ) 

REAC (5,67) ( A( I ),E0( I ) ,GAR!N( I),GAMG( I ),C-ARE( I ) ,1 = 1 ,NRES) 
REAC (5,71) NUM8G 
REAC (5,53) ( GG ( I ) , I* 1, NUR1BG ) 

REAC(5,75 ) 12,14,16,18 
WRITE (6,54) 

WRITE (6,50 (TITLE(R),M = 1,24) 



WRITE (6,55) 

WRITE (6,56) T, ANU,UKT,SIGCT,SIGFT,SIGMAT,SIGPOT 
DC 2 K= 1 , NCU , 52 
IHI=NINO (NCU,K + 51 ) 

C NEW PAGE 

WRITE (6,54) 

WRITE (6,57) 

WRITE (6,58) (J,l!G(J),DU(J),J=K,IHI) 

2 CONTINUE 

DC 3 1=1, N RES 

GAM ( I )=GAMN (I )+GAKG< I )+GAME( I ) 

CCN'1= ( 2.607 E06/(G AM ( I )**2) )*G*GAMN( I ) 

CCN2 = ( A ( I ) +1 .C ) / A ( I ) 

CCN3 ( I )=1.C/C0N2 
AEC = AE S ( EC ( I ) ) 

SIGCCII )= (CCN1*GAMG( l)*C0N2 **1.5 > /SQRT ( AEO ) 

S I GFC ( I ) * ( G AME( I ) *S IGCO ( I) )/GAMG{ I ) 

SIGSC ( I )= (CCN1*GAMN( I )*C0N2 **2)/AE0 

SIC-TC (I ) = 2.607E06*G*GAMN( I ) / ( AEQ *G AM ( I ) > 

3 CCNT IMJE 

DC 4 1= 1 *NRES ,52 
IH I = N INO (NPES ,1+51 ) 

C NEW PAGE 

WRITE (6,54) 

WRITE (6,55) 

WRITE (6,60) (J,A(J),EO(J),GAMN(J),GAMG(J),GAME(J),J=I,IHI) 

4 CONTINUE 

DC 5 1=1, N RES, 52 
IHI=NINC (NRES,I+51) 

C NEW PAGE 

WRITE (6,54) 

WRITE (6,61) 

WRITE (6,62) (J»EQ(J)»SI GCO (J),SIGFO(J) , S I GSO ( J ) ,SIGTO(J) ,J 

5 CONTINUE 
WRITE (6,55) 

C 

NCCUNT=0 
NT AL L V= 1 
KKK = C 
JJJ = C 
C 

ET=EXP(-LKT)*1.0E07 

6 CC 13 1=1, N RES 
IF ( KKK ) 7,7,10 

7 X ( I ) = (2 .C * ( ET-EO ( I ) ) *CON 3 ( I ) ) /G AM ( I ) 

IF (T) 8,8,5 

8 TCCF=C.O 
GC TC 11 

9 TCCF=293 .0 
GC TC 11 

10 TCCF=T 

11 THETA(I) = GAM( I )/C0N3( I )*SGRT( A( I )/(TDOP*ABS(EO( I ) )*3.446E 

13 CCNT INU E 

14 CALL PSIFHI 

IF (KKK) 25,15,25 

15 SIGCTP=C.O 
SIGFTP=0.0 
ZEC1=SQRT (ET) 

CC 16 1 = 1 ,NRES 


I , IHI ) 


04 ) ) 
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SIGCTP=SIGCTP+SIGCO( I )*PSI ( I ) 
SIGFTP=SIGFTP+SIGFG( I )*PSI ( I ) 

16 CONTINUE 
SCV=SIGC7-S IGCTP/ZE01 
SFV = S IGFT-S IGFT P/ZED1 
AKC = SCV *2 EC 1 

AK F= S FV*Z EC 1 

WRITE (6,64) SCV,SFV 

A E N = C • 0 

DC 17 1=1 , N R ES 
AEN = AEN + MI) 

17 CONTINUE 

SIGP=.27154671*( AEN/FLOAT(NRES) )**. 66666667 
IF (SCV.LT.O.O.ANC.SFV.LT.O.C) GO TO 1 
I F ( 1 2 ) 19,19,18 

18 SIGP=SIGFCT 
GC TC 20 

19 CONTINUE 

WRITE (6,65) SIGP 
IF (SIGP) 1,20,20 
2C WR I T E ( 6,54) 

WR IT E ( 6,66) 

LA=1 

I F ( 1 6 ) 22,22 ,21 

21 E=CSTART 
U=16.118C96-AL0G(E) 

GC TC 23 

22 U=LST ART 
E=1.CE07*EXP(-U) 

23 DC 24 1=1 » N R E S 

24 X ( I ) = (2.0*CON3( I)/GAM( I) )*(E-EO( I) ) 
KKK=KKK+1 
GC TC 6 
C 

25 SIGCN=O.C 
SIGFN=O.C 
SIC-SN=O.C 
ZED 1 = SQRT ( E ) 

DC 26 1=1 ,NRES 
ZILCH1=SCRT(SIGSC( I ) ) 

Z I LCF2= PF I ( I )*SGRT(G*SIGP*SIGSO( I ) ) 

ZILCF3 = PS Id) 

SIG1 ( I ) =FH I ( I )*Z ILCH1 
SIG2 ( I ) = 2 I LCH2 
SIG2(I)=ZILCH3*ZILCH1 
WCRK = Z I LCH2/CCN3 ( I ) 

TEEP=SIGSC( I ) *PS I ( I ) 

SIGCN = S ICCN + S IGCC ( I >*Z I LCH3 
SIGFN=SICFN+SIGFC( I )*ZILCH3 
SIGSN=SIGSN+TEMP+WCRK 

26 CONTINUE 

SIGCN=( SIGCN + AKC )/ZECl 
SIGFN=( SIGFN + AKF ) /ZEC1 
SIGSN=SIGSN+SIGP 
I F ( 1 8 ) 30,30,27 

27 SIGSN=S IGP 


29 



28 

29 

30 

31 

32 

C 


33 


34 


C 

37 

38 

40 


41 


42 


C 

44 

45 


43 
4 6 


47 


46 


DC 29 1 = 1 tMJWBG 
SIGY1=0. 

SIGY2=0. 

S1GY3=0. 

CC 28 J=1,NRES 

IF (A8S(G-GG( I) ) .GT.C.00C01) GO TO 28 
SIGY1=SICY1*SIG1(J) 

SIGY2 = SIC-Y2 + SIG2(J) 

SIGY3=SIGY3+SIG3(J) 

CONTINUE 

SIGSN=SIGSN+0.25*SIGY1**2+SIGY2+SIGY3**2 

CONTINUE 

IF (SIGSM 31,32,32 
JJJ=1 

$IGA=SIGCN+S IGFN 
ANLSIG=ANU*SIGFN 

IF (NCOUM-54) 34,34,33 
WRITE! 6,54) 

WRITE! 6,66) 

NCCLN T=G 

WRITE! 6,68) NT A L L Y , U , E , S I G A , S I GSN , ANUS I G 

IWCPC = NWCRCS-NT ALLY + 1 

NTALLY=N1 ALLY+1 

NCCLNT=NCCLNT+1 

DATLNS ( IWCPC , 1 ) = E 

OATLNS! IWCRC, 2) = SICA 

DATLNS ! IWCRC., 3) = SIGSN 

I F ! 1 6 ) 44,44,37 

IF ! ! UG C L A )-E) + O.OOOCl) 38,40,40 
E=E-CU( LA ) 

U=16.118C96-ALOG(E) 

GC TC 43 
E=UG (LA) 

LA=L A+l 

IF ( (LG(NCU)-E)+.OOOC01 ) 38,41,41 

CONTINUE 

IF(JOJ) 1,1,42 

WRITE (6,54) 

WRITE (6,70) 

WRITE(18,7C) 

GC TC 1 

IF ( (LG(LA)-U)-O.COOCl) 47,47,45 
U=L+CU ( LA ) 

E=1.CE07*EXP(-U ) 

CONTINUE 
CC 46 1=1 ,NRES 

X! I ) = (2.C4C0N3! I ) /G AM ( I) ) * ( E-EO ( I )) 

GC TC 14 
U=LG ( L A ) 

LA = L A *1 

IF ( (UG(NCU)-U)-.COOGl) 48,46,45 

CONTINUE 

IF (JJJ) 1,1,49 


30 



49 WRITE (6,54) 

WRITE (6,70 
GC TC 1 

C 

50 FORMAT ( 1 2 £6 ) 

51 FCRMT ((FIG. 5, 215) 

52 FCRMT (JF10.6) 

53 FCRMT (7F1C.7) 

54 FCRMT ( 1H 1 ) 

55 FCRMT ( 1 HK ) 

56 FCRMT ( 3 X , 6HTEMP = F 1 2. 2 , 3X , 4HNU =F6 . 3 , 3X , (HU KT =F 6 . 4 , 3 X , 1 2HS T G C 

1AP KT =F10.4,3X, 10HSIG F KT =F1C.4,3X,9HSIG TOT = F 1 C . 4/3 X , 9H SI G PO 
2T = F9 .3 ) 

57 FCRMT ( 10X 1 2H K,7X,3H UG,11X,3H OU ) 

58 FCRMT (112, 2 FI 3. 5) 

59 FCRMT ( 1HJ/3X, 1HI , 15H ATOMIC WEIGHT, 3X, 1C HR ES ENERGY, 6X ,7HGAMMA 
1N,7X,7HGAMMA G, 7X , 7FG AMM A F) 

60 FORMAT ( 14 , F 14. 3 , F 18 . 4, 3F14. E ) 

61 FORMAT (4X»1HI»3X,10HRES ENERGY , 5X , 5HS IGCO , 9X , 5HS IGFO , 9X , 5HS I GSO , 9 
IX, 5HSIGTC ) 

62 FORMAT ( 1 5 , FI 3. 4 ,4F14 .4 ) 

64 FORMAT (10X,10HSIC C /V = , F 1 2 . 4 , 5X , 1CHS I G F /V =,F12.4) 

65 FORMAT (1H /10X,10HSIG POT =F12.4) 

66 FORMAT ( 6 X , 1HK, 6X , 1HU , 1 1 X , 1H E , 1 1 X , 5HS I G A,8X,6HSIG ES,7X,8HNU SIG 
IF) 

67 FORMAT ( 1 OX , F 10 . 7 , 10X , 4F 10 . 7 ) 

68 FORMAT ( 1 7 , 2X ,F8 . 5 , F 1 3 . 5 , 6F 1 3 .4 ) 

70 FORMAT (10X.22HSIG S IS NEG AT SOME E) 

71 FORMAT (15) 

72 FORM AT ( IX, II, IX, II, IX, II, IX, ID 
ENC 
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